# data handling
library(readxl)
library(magrittr)
library(reshape2)
library(plyr)
library(rstatix)
##
## Attaching package: 'rstatix'
## The following objects are masked from 'package:plyr':
##
## desc, mutate
## The following object is masked from 'package:stats':
##
## filter
# general R plots
library(ggplot2)
library(ggpubr)
##
## Attaching package: 'ggpubr'
## The following object is masked from 'package:plyr':
##
## mutate
library(pheatmap)
# survival plots
library(survival)
library(survminer)
##
## Attaching package: 'survminer'
## The following object is masked from 'package:survival':
##
## myeloma
# figure outdir
outdir <- "~/Figures/"
# colonisation
col_palette <- c("Clean"="darkblue", "Non_col"="grey")
# isolate colours
iso_select_1 <- c("PBS",
"Enterocloster clostridioformis",
"Escherichia coli str. A7",
"Escherichia coli str. Fergi",
"Anaerostipes sp000508985",
"S. Typhimurium")
iso_select_2 <- c("E. clostridioformis",
"E. coli A7",
"E. coli Fergi",
"A. faecis")
iso_pal_1 <- get_palette("Set1", 6)
iso_pal_2 <- get_palette("Set1", 6)[c(2:5)]
names(iso_pal_1) <- iso_select_1
names(iso_pal_2) <- iso_select_2
iso_pal <- c(iso_pal_1, iso_pal_2, "LB"="#E41A1C")
phycol=c("#377eb8","#8CB302","#008C74","#d95f02","#FF974F","#FFED6F","#FDCDAC","#ffd92f","#e22426","#B3B3B3","#FBB4AE","#984ea3","#35478C","#7FC97F","#FF73C5","#BF5693")
phylabs=c("Actinobacteriota","Bacteroidota","Campylobacterota","Cyanobacteria","Deferribacterota","Desulfobacterota","Elusimicrobiota","Firmicutes","Firmicutes_A","Firmicutes_B","Firmicutes_C","Proteobacteria","Spirochaetota","Thermotogota","Verrucomicrobiota","Verrucomicrobiota_A")
names(phycol) <- phylabs
classcol=c("#377eb8","#8CB302","#008C74","#d95f02","#FF974F","#FFED6F","#FDCDAC","#ffd92f","#e22426","#B3B3B3","#FBB4AE","#984ea3","#35478C","#7FC97F","#FF73C5","#BF5693")
classlabs=c("Actinomycetia","Bacteroidia","Campylobacterota","Cyanobacteria","Deferribacterota","Desulfobacterota","Elusimicrobiota","Bacilli","Clostridia","Firmicutes_B","Firmicutes_C","Gammaproteobacteria","Spirochaetota","Thermotogota","Verrucomicrobiota","Verrucomicrobiota_A")
names(classcol) <- classlabs
dp = "Figure data.xlsx"
readxl::excel_sheets(dp)
## [1] "GF_screen_data" "Post_STm_oedema"
## [3] "Post_STm_histology_scores" "in_vivo_CFU"
## [5] "Colonisation_resistance" "Spent_media_glucose"
## [7] "Pre_STm_mucus_scores" "IEC_RTqPCR_data"
## [9] "GO_cell_cycle" "GO_leukocyte_activation"
## [11] "Flow_data_IEL" "Flow_data_nonIEL"
## [13] "SPF_screen_data"
GF screen of MCC isolates
# read data
d <- read_excel(path = dp, sheet = "GF_screen_data")
# fix taxonomy
# define species
d$Species <- unlist(lapply(strsplit(d$Commensal_taxonomy_r95, split = ";.__"),
function(x) {
x[7]
}))
d$Species[is.na(d$Species)] <- "PBS"
# define class
d$Class <- unlist(lapply(strsplit(d$Commensal_taxonomy_r95, split = ";.__"),
function(x) {
x[3]
}))
d$Class[is.na(d$Class)] <- "PBS"
# define phylum
d$Phylum <- unlist(lapply(strsplit(d$Commensal_taxonomy_r95, split = ";.__"),
function(x) {
x[2]
}))
d$Phylum[is.na(d$Phylum)] <- "PBS"
# clean up data
d[,grep(pattern = "Weight_", colnames(d))] <- apply(X = d[,grep(pattern = "Weight_", colnames(d))],
MARGIN = 2,
FUN = function(x) {
x[x == 0] <- NA
as.numeric(x)
}) %>% as.data.frame
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
# define species print names
d$Species_id <- d$Species
d$Species_id[grep(pattern = "coli", d$Species)] <- paste(d$Species[grep(pattern = "coli", d$Species)], "str.", d$Treatment[grep(pattern = "coli", d$Species)])
# generate percentage weight loss
d_perc <- apply(X = d[,grep(pattern = "Weight_", colnames(d))],
MARGIN = 1,
function(x) {
x/x[1]*100
}) %>% t %>% as.data.frame
colnames(d_perc) <- paste0(colnames(d_perc), "_perc")
d <- cbind(d, d_perc)
colnames(d)[grep(pattern = "_perc", colnames(d))] <- c("dpi0", "dpi1", "dpi2", "dpi3", "dpi4")
# edit palette
classpal <- classcol[names(classcol) %in% unique(d$Class)]
## plot 1dpi weights
i <- unlist(lapply(split(d$dpi1, f = d$Species_id), median))
dpi1_order <- names(i)[order(i, decreasing = TRUE)]
p <- ggboxplot(d, x = "Species_id", y = "dpi1",
order = rev(dpi1_order),
outlier.size = 0.3,
fill = "Class", palette = classpal,
ylab = FALSE) +
ylim(c(80,110))
p1 <- ggpar(p, rotate = TRUE, legend = "none", ylab = "weight at 1 dpi (%)") +
stat_compare_means(ref.group = "PBS",
label = "p.signif",
hide.ns = TRUE,
size = 3, vjust = 0.8)
p1 <- p1 + theme(axis.ticks.y = element_blank(),
axis.line.y = element_blank(),
axis.text.y = element_blank()) +
theme(axis.title.x = element_text(vjust=19))
## plot survival
# get data
d$Survival_day <- apply(X=d[,grep(pattern = "^dpi", colnames(d))], MARGIN = 1, function(x) {
colnames(d)[grep(pattern = "^dpi", colnames(d))][length(x) + 1 - as.numeric(which(!is.na((rev(x))))[1])]
}) %>% gsub(pattern = "dpi", replacement = "") %>% as.numeric
# summarise the data -> mean_se
ds <- ddply(.data = d, .variables = c("Species_id", "Class"), .drop = TRUE,
.fun = function(xx, col) {
c(mean = mean(xx[[col]]),
se = sd(xx[[col]]) / sqrt(length(xx[[col]]))
)
},
"Survival_day"
)
d2 <- merge(x = d, y = ds, by = c("Species_id", "Class"), all.x = TRUE)
p <- ggline(d2, x = "Species_id", y = "Survival_day", color = "Class",
palette = classpal,
group = "Species_id",
plot_type = "p",
add = "mean",
#size = 0.4,
order = rev(dpi1_order),
ylab = FALSE) +
geom_errorbar(aes(ymin=mean-se, ymax=mean+se), size = 0.4, width = 0.1, colour="black") +
geom_point(aes(y = mean, color = Class), size = 2) +
ylim(c(1,4)) +
stat_compare_means(ref.group = "PBS", label = "p.signif", hide.ns = TRUE,
size = 3, vjust = 0.8)
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
p2 <- ggpar(p, rotate = TRUE, legend = "right",
legend.title = "Taxonomic class",
ylab = "survival (dpi)")
p2 <- p2 +
theme(axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
axis.line.y = element_blank()) +
theme(axis.title.x = element_text(vjust=19))
## plot colonisation status
d$colstat <- 1
# change PBS to being non-colonised
d$Col_status[d$Treatment == "PBS"] <- "Non_col"
pb <- ggbarplot(d, x = "Species_id", y = "colstat",
add = "mean",
order = rev(dpi1_order),
fill = "Col_status",
palette = c("black", "white"),
ylab = FALSE,
xlab = FALSE)
pb <- ggpar(pb, rotate = TRUE, legend = "none", x.text.angle = 30) +
scale_y_continuous(breaks=c(0.5), labels = "Colonisation status") # single central tick mark
pb <- pb + theme(axis.ticks.y = element_blank(),
axis.line = element_blank())
ggarrange(pb, p1, p2, nrow = 1, widths = c(0.252,0.3,0.4), align = "h")
## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_boxplot()`).
## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_compare_means()`).
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's fill values.
# weight loss at 1dpi
wilcox_test(data = d, formula = dpi1 ~ Treatment, ref.group = "PBS", p.adjust.method = "BH")
## # A tibble: 18 × 9
## .y. group1 group2 n1 n2 statistic p p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <chr>
## 1 dpi1 PBS A12 18 4 31 0.712 7.84e-1 ns
## 2 dpi1 PBS A14 18 7 45.5 0.303 6.06e-1 ns
## 3 dpi1 PBS A17 18 3 9 0.08 2.85e-1 ns
## 4 dpi1 PBS A2 18 2 5 0.126 3.24e-1 ns
## 5 dpi1 PBS A22 18 3 31 0.74 7.84e-1 ns
## 6 dpi1 PBS A26 18 5 33 0.403 6.18e-1 ns
## 7 dpi1 PBS A43 18 5 34 0.446 6.18e-1 ns
## 8 dpi1 PBS A61 18 8 89 0.367 6.18e-1 ns
## 9 dpi1 PBS A68 18 4 44 0.538 6.92e-1 ns
## 10 dpi1 PBS A7 18 3 0 0.002 9 e-3 **
## 11 dpi1 PBS A82 18 9 40 0.035 1.59e-1 ns
## 12 dpi1 PBS A92 18 18 0 0.00000000022 3.96e-9 ****
## 13 dpi1 PBS B37 18 3 25 0.887 8.87e-1 ns
## 14 dpi1 PBS C48 18 4 29.5 0.609 7.31e-1 ns
## 15 dpi1 PBS C49 18 9 48 0.095 2.85e-1 ns
## 16 dpi1 PBS C57 18 4 21 0.227 5.11e-1 ns
## 17 dpi1 PBS Fergi 18 4 0 0.000273 2 e-3 **
## 18 dpi1 PBS VP_Ef 18 4 26 0.434 6.18e-1 ns
# survival
wilcox_test(data = d, formula = Survival_day ~ Treatment, ref.group = "PBS", p.adjust.method = "BH")
## # A tibble: 18 × 9
## .y. group1 group2 n1 n2 statistic p p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <chr>
## 1 Survival_day PBS A12 18 4 44 3.4 e-1 5.56e-1 ns
## 2 Survival_day PBS A14 18 7 26 9 e-3 3.4 e-2 *
## 3 Survival_day PBS A17 18 3 33 4.17e-1 5.56e-1 ns
## 4 Survival_day PBS A2 18 2 13 4.5 e-1 5.56e-1 ns
## 5 Survival_day PBS A22 18 3 24 7.34e-1 8.26e-1 ns
## 6 Survival_day PBS A26 18 5 37 4.63e-1 5.56e-1 ns
## 7 Survival_day PBS A43 18 5 28 1.23e-1 3.69e-1 ns
## 8 Survival_day PBS A61 18 8 61 4.48e-1 5.56e-1 ns
## 9 Survival_day PBS A68 18 4 35 9.53e-1 1 e+0 ns
## 10 Survival_day PBS A7 18 3 0 1 e-3 1.3 e-2 *
## 11 Survival_day PBS A82 18 9 81 1 e+0 1 e+0 ns
## 12 Survival_day PBS A92 18 18 4 1.9 e-7 3.42e-6 ****
## 13 Survival_day PBS B37 18 3 33 4.17e-1 5.56e-1 ns
## 14 Survival_day PBS C48 18 4 26 2.95e-1 5.56e-1 ns
## 15 Survival_day PBS C49 18 9 99 1.44e-1 3.7 e-1 ns
## 16 Survival_day PBS C57 18 4 26 2.95e-1 5.56e-1 ns
## 17 Survival_day PBS Fergi 18 4 6 3 e-3 1.3 e-2 *
## 18 Survival_day PBS VP_Ef 18 4 6 3 e-3 1.3 e-2 *
# subset data
iso_select <- c("PBS","Enterocloster clostridioformis")
df <- d[d$Species_id %in% iso_select,]
df$Species_id <- factor(df$Species_id, levels = iso_select)
# format data
df$idvar <- paste0("M", seq(1, nrow(df), 1))
dfm <- melt(df, measure.vars = colnames(df)[grep(pattern = "^dpi", colnames(df))])
dfm$dpi <- gsub(pattern = "dpi", replacement = "", dfm$variable) %>% as.numeric
# weight loss dpi
ddply(.data = df, .variables = "Species_id",
.fun = function(x) {
c(N=nrow(x),
apply(x[,grep("^dpi", colnames(x))],
MARGIN = 2,
FUN = function(y) {
sum(!is.na(y))
}))
})
## Species_id N dpi0 dpi1 dpi2 dpi3 dpi4
## 1 PBS 18 18 18 4 0 0
## 2 Enterocloster clostridioformis 18 18 18 18 16 8
# remove timepoints where mice have died -> plot only up to day 3
dfm_v2 <- rbind(dfm[dfm$dpi %in% c(0:1) & dfm$Treatment == "PBS",],
dfm[dfm$dpi %in% c(0:3) & dfm$Treatment == "A92",])
p <- ggline(dfm_v2, x = "dpi", y = "value",
group = "Species_id", color = "Species_id", plot_type = "b",
add = "mean_se",
size = 1,
palette = "Set1",
xlab = "days post-infection",
ylab = "weight (% of day 0)") +
geom_hline(yintercept = 80, linetype = "dashed") +
theme(legend.title = element_blank()) +
stat_compare_means(label = "p.format", label.x = 3.5, label.y = 97) +
stat_compare_means(label = "p.signif", label.x = 2.9, label.y = 97)
p_dpi2 <- ggpar(p, legend = "right")
p_dpi2
## Warning: Removed 2 rows containing non-finite outside the scale range
## (`stat_summary()`).
## Warning: Removed 2 rows containing non-finite outside the scale range
## (`stat_compare_means()`).
## Removed 2 rows containing non-finite outside the scale range
## (`stat_compare_means()`).
## survival
df$Status <- 2
surv_diff <- survdiff(Surv(Survival_day, Status) ~ Species_id, data = df)
d_survival <- ddply(.data = df, .variables = "Species_id",
.fun = function(x) {
dpi0_surv <- length(which(!is.na(x$dpi1)))
dpi1_surv <- length(which(!is.na(x$dpi2)))
dpi2_surv <- length(which(!is.na(x$dpi3)))
dpi3_surv <- length(which(!is.na(x$dpi4)))
dpi0 <- dpi0_surv/dpi0_surv*100
dpi1 <- dpi1_surv/dpi0_surv*100
dpi2 <- dpi2_surv/dpi0_surv*100
dpi3 <- dpi3_surv/dpi0_surv*100
c(dpi0=dpi0,
dpi1=dpi1,
dpi2=dpi2,
dpi3=dpi3,
dpi4=0)
}
)
d_survival[,grep("dpi", colnames(d_survival))] <- apply(d_survival[,grep("dpi", colnames(d_survival))],
MARGIN = 1,
FUN = function(x) {
i <- sum(x == 0) - 1
if ( i != 0) {
x[6-c(1:i)] <- NA
}
as.numeric(x)
}) %>% t
dm_survival <- melt(data = d_survival, id.vars = "Species_id")
dm_survival$dpi <- gsub(pattern = "dpi", replacement = "", dm_survival$variable) %>%
as.numeric
p <- ggplot(data = dm_survival, mapping = aes(x = dpi, y = value, group = Species_id)) +
geom_step(aes(color = Species_id), size = 1) +
geom_point() +
theme_pubr() +
theme(legend.title = element_blank())
psurv <- ggpar(p, palette = "Set1", xlab = "days post-infection", ylab = "percent survival",
legend = "right") +
annotate(geom="text", label = "**** p = 5e-09", x=3.3, y=98)
psurv
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_step()`).
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).
Histology
oedema <- read_excel(path = dp, sheet = "Post_STm_oedema")
colnames(oedema)[3] <- "oedema_perc"
colnames(oedema)[10] <- "Sample"
oedema$oedema_perc <- as.numeric(oedema$oedema_perc)
## Warning: NAs introduced by coercion
oedema_score_SS <- ddply(oedema[oedema$Scorer == "SS",], .variables = c("Tissue", "Sample", "Group"),
.fun = function(x) {
m <- mean(x$oedema_perc)
if (m >= 80) {
s <- 3
} else if (m >= 50) {
s <- 2
} else if (m >= 20) {
s <- 1
} else {
s <- 0
}
c(Mean=m, Oedema_score=s)
})
oedema_score_KR <- ddply(oedema[oedema$Scorer == "KR",], .variables = c("Tissue", "Sample", "Group"),
.fun = function(x) {
m <- NA
s <- median(x$`submucosal oedema`)
c(Mean=m, Oedema_score=s)
})
oedema_score <- ddply(rbind(oedema_score_SS, oedema_score_KR), .variables = c("Tissue", "Sample", "Group"),
.fun = function(x) {
m <- mean(x$Oedema_score)
c(Oedema_score=m)
})
# read histology data
histo <- read_excel(path = dp, sheet = "Post_STm_histology_scores")
histo$Group[histo$Group == "A92"] <- "E. clostridioformis"
# format data
histo$`Neuts/mm2` <- as.numeric(histo$`Neuts/mm2`)
## Warning: NAs introduced by coercion
histo$mononuclear_infiltrate <- as.numeric(histo$mononuclear_infiltrate)
histo$epithelial_integrity <- as.numeric(histo$epithelial_integrity)
histo$cryptitis <- as.numeric(histo$cryptitis)
histo$goblet_cells_per_section <- as.numeric(histo$goblet_cells_per_section)
## Warning: NAs introduced by coercion
histo$neutrophils_per_section <- as.numeric(histo$neutrophils_per_section)
## Warning: NAs introduced by coercion
# HPF area = 262.1 um x 262.1 um
HPF_area_mm2 = 0.06869641
# calculate average scores for samples
histo_scores <- ddply(.data = histo, .variables = c("Sample_index", "Tissue", "Sample", "Group"),
.fun = function(x) {
Mean = mean(x$`Neuts/mm2`, na.rm = TRUE)
c(Neuts_count = (Mean * HPF_area_mm2) / 2, # normalise for HPF area
GC_count = mean(as.numeric(x$goblet_cells_per_section), na.rm = TRUE),
Mononuclear_infiltrate=median(as.numeric(x$mononuclear_infiltrate), na.rm = TRUE),
Epithelial_integrity=median(as.numeric(x$epithelial_integrity), na.rm = TRUE),
Cryptitis=median(as.numeric(x$cryptitis), na.rm = TRUE)
)
})
# Neutrophils
histo_scores$Neutrophil_score <- sapply(histo_scores$Neuts_count, function(x) {
x <- as.numeric(x)
if (is.na(x)) {
s=NA
} else if (x > 100) {
s=4
} else if (x > 60) {
s=3
} else if (x > 20) {
s=2
} else if (x >= 5) {
s=1
} else {
s=0
}
return(s)
})
# Goblet cells
histo_scores$GC_score <- sapply(histo_scores$GC_count, function(x) {
x <- as.numeric(x)
if (is.na(x)) {
s=NA
} else if (x > 28) {
s=0
} else if (x > 10) {
s=1
} else if (x >= 1) {
s=2
} else if (x < 1) {
s=3
}
return(s)
})
# generate final scores
histo_scores_final <- ddply(histo_scores, .variables = c("Sample", "Tissue", "Group"),
.fun = function(x) {
c(Mononuclear_infiltration = mean(x$Mononuclear_infiltrate, na.rm = TRUE),
Epithelium = mean(x$Epithelial_integrity, na.rm = TRUE),
Cryptitis = mean(x$Cryptitis, na.rm = TRUE),
PMNs = mean(x$Neutrophil_score, na.rm = TRUE),
GCs = mean(x$GC_score, na.rm = TRUE)
)
})
histo_m <- melt(data = histo_scores_final,
id.vars = c("Sample", "Tissue", "Group"),
measure.vars = c("Mononuclear_infiltration",
"Cryptitis",
"Epithelium",
"PMNs",
"GCs"))
# combine with oedema data
oedema_m <- melt(oedema_score, measure.vars = "Oedema_score")
histo_scores_final <- rbind(oedema_m, histo_m)
histo_scores_final$variable <- droplevels(histo_scores_final$variable) %>% as.character
histo_scores_final$variable[histo_scores_final$variable == "Oedema_score"] <- "Oedema"
# set levels
histo_scores_final$Group <- factor(histo_scores_final$Group, levels = unique(histo_scores_final$Group))
# define names
histo_scores_final$Obs <- NA
histo_scores_final$Obs[histo_scores_final$variable == "Oedema"] <- "Submucosal oedema"
histo_scores_final$Obs[histo_scores_final$variable == "Mononuclear_infiltration"] <- "Mononuclear infiltrate"
histo_scores_final$Obs[histo_scores_final$variable == "Epithelium"] <- "Epithelial integrity"
histo_scores_final$Obs[histo_scores_final$variable == "Cryptitis"] <- "Cryptitis"
histo_scores_final$Obs[histo_scores_final$variable == "PMNs"] <- "PMN infiltration"
histo_scores_final$Obs[histo_scores_final$variable == "GCs"] <- "Goblet cells"
histo_total <- ddply(histo_scores_final, .variables = c("Sample", "Group", "Tissue"),
.fun = function(x) {
c(Total_score=sum(x$value, na.rm = TRUE))
})
p_hist <- ggbarplot(histo_total, x = "Tissue", y = "Total_score",
fill = "Group", palette = "Set1",
order = c("Ileum", "Caecum", "Colon"),
ylab ="Histology score",
xlab = FALSE,
add = c("mean_se"), position = position_dodge(0.8))
p_hist <- ggpar(p_hist, legend = "right") +
theme(legend.title = element_blank()) +
scale_y_continuous()
p_hist
p <- ggbarplot(histo_scores_final[histo_scores_final$Tissue == "Caecum",], x = "Group", y = "value",
fill = "Obs", palette = "Greys",
add = c("mean_se"),
xlab = FALSE, ylab = "Score")
p_hcaecum <- ggpar(p, legend = "right") +
theme(legend.title = element_blank(),
axis.title.x = element_blank())
p_hcaecum
In vivo S.Tm titres.
# get data
vivo_cfu <- read_excel(path = dp, sheet = "in_vivo_CFU")
# define species
vivo_cfu$Species <- unlist(lapply(strsplit(vivo_cfu$Commensal_taxonomy_r95, split = ";.__"),
function(x) {
x[7]
}))
vivo_cfu$Species[is.na(vivo_cfu$Species)] <- "PBS"
vivo_cfu$Species[vivo_cfu$Species == "Enterocloster clostridioformis"] <- "E. clostridioformis"
# change names for publication
vivo_cfu$Group <- vivo_cfu$Species
vivo_cfu$Group[vivo_cfu$Group == "PBS"] <- "GF"
vivo_cfu$Group[vivo_cfu$Group == "Anaerostipes sp000508985"] <- "A. faecis"
vivo_cfu$Group[vivo_cfu$Group == "Escherichia coli"] <- "E. coli"
# add levels
vivo_cfu$Group <- factor(vivo_cfu$Group, levels = unique(vivo_cfu$Group))
# create palette
monocol_pal <- get_palette("Set1", 4)
names(monocol_pal) <- unique(vivo_cfu$Group)
# default comparisons for stats
monocol_comp = list(c("GF", "E. clostridioformis"), c("GF", "A. faecis"), c("GF", "E. coli"))
p_caecal_contents <- ggboxplot(vivo_cfu, x = "Group", y="Caecal_content_STm_CFUg",
fill = "Group",
palette = monocol_pal,
order = unique(vivo_cfu$Group),
size = 1,
ylab = "STm CFU/g caecal contents",
xlab = "Caecal contents",
add = c("dotplot")) +
stat_compare_means(comparisons = monocol_comp,
label = "p.signif",
bracket.size = 0.5,
tip.length = 0)
p_caecal_contents <- ggpar(p_caecal_contents, legend = "right") +
scale_y_log10() +
theme(axis.ticks.x = element_blank(),
axis.text.x = element_blank(),
legend.title = element_blank())
p_caecal_contents
## Bin width defaults to 1/30 of the range of the data. Pick better value with
## `binwidth`.
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
# correct stats for multiple copmarisons
wilcox_test(vivo_cfu, formula = Caecal_content_STm_CFUg ~ Group, p.adjust.method = "BH", ref.group = "GF")
## # A tibble: 3 × 9
## .y. group1 group2 n1 n2 statistic p p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <chr>
## 1 Caecal_conte… GF E. cl… 8 13 85 1.6 e-2 2.4 e-2 *
## 2 Caecal_conte… GF A. fa… 8 7 35 4.63e-1 4.63e-1 ns
## 3 Caecal_conte… GF E. co… 8 7 56 3.11e-4 9.33e-4 ***
# get fold change values for E. coli vs E. clostridioformis
caecum_medians <- ddply(.data = vivo_cfu, .variables = c("Group"), .fun = function(x) {
median(x$Caecal_content_STm_CFUg)
})
caecum_medians
## Group V1
## 1 GF 9.25e+08
## 2 E. clostridioformis 3.16e+08
## 3 A. faecis 5.97e+08
## 4 E. coli 2.40e+07
# E. coli
caecum_medians$V1[1]/caecum_medians$V1[4]
## [1] 38.54167
# 38.54167
# E. clostridioformis
caecum_medians$V1[1]/caecum_medians$V1[2]
## [1] 2.927215
# 2.927215
vivo_cfu$mLN_STm_CFUg[vivo_cfu$mLN_STm_CFUg == 0] <- 1
p_mLN <- ggboxplot(vivo_cfu, x = "Group", y="mLN_STm_CFUg",
fill = "Group",
palette = monocol_pal,
size = 1,
ylab = "STm CFU/g mLN",
xlab = "mLN",
add = c("dotplot")) +
stat_compare_means(comparisons = monocol_comp,
label = "p.signif",
# label.y = 5.9,
# size = 5,
bracket.size = 0.5,
tip.length = 0)
p_mLN <- ggpar(p_mLN, legend = "right") +
yscale("log10") +
scale_y_log10() +
theme(axis.ticks.x = element_blank(),
axis.text.x = element_blank(),
legend.title = element_blank())
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
p_mLN
## Bin width defaults to 1/30 of the range of the data. Pick better value with
## `binwidth`.
## Warning in wilcox.test.default(c(1.07918124604762, 0, 0, 5.20411998265593, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(c(1.07918124604762, 0, 0, 5.20411998265593, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(c(1.07918124604762, 0, 0, 5.20411998265593, :
## cannot compute exact p-value with ties
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
# correct stats for multiple copmarisons
wilcox_test(vivo_cfu, formula = mLN_STm_CFUg ~ Group, p.adjust.method = "BH", ref.group = "GF")
## # A tibble: 3 × 9
## .y. group1 group2 n1 n2 statistic p p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <chr>
## 1 mLN_STm_CFUg GF E. clostri… 8 13 70 0.197 0.296 ns
## 2 mLN_STm_CFUg GF A. faecis 8 7 27 0.954 0.954 ns
## 3 mLN_STm_CFUg GF E. coli 8 7 49 0.007 0.022 *
# make data numerics and fix 0's for log range
vivo_cfu$Liver_STm_CFUg <- as.numeric(vivo_cfu$Liver_STm_CFUg)
## Warning: NAs introduced by coercion
vivo_cfu$Liver_STm_CFUg[vivo_cfu$Liver_STm_CFUg == 0] <- 1
p_liver <- ggboxplot(vivo_cfu, x = "Group", y="Liver_STm_CFUg",
fill = "Group",
palette = monocol_pal,
size = 1,
ylab = "STm CFU/g liver",
xlab = "Liver",
add = c("dotplot")) +
stat_compare_means(comparisons = monocol_comp,
label = "p.signif",
#label.y = 7.6,
#size = 5,
bracket.size = 0.5,
tip.length = 0)
p_liver <- ggpar(p_liver, legend = "right") +
scale_y_log10() +
theme(axis.ticks.x = element_blank(),
axis.text.x = element_blank(),
legend.title = element_blank())
p_liver
## Warning: Removed 6 rows containing non-finite outside the scale range
## (`stat_boxplot()`).
## Bin width defaults to 1/30 of the range of the data. Pick better value with
## `binwidth`.
## Warning: Removed 6 rows containing missing values or values outside the scale range
## (`stat_bindot()`).
## Warning: Removed 6 rows containing non-finite outside the scale range
## (`stat_signif()`).
## Warning in wilcox.test.default(c(4.22788670461367, 1.2380461031288, 0,
## 2.46389298898591, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(c(4.22788670461367, 1.2380461031288, 0,
## 2.46389298898591, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(c(4.22788670461367, 1.2380461031288, 0,
## 2.46389298898591, : cannot compute exact p-value with ties
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
# correct for multiple comparisons
wilcox_test(vivo_cfu, formula = Liver_STm_CFUg ~ Group, comparisons = monocol_comp, p.adjust.method = "BH", ref.group = "GF")
## # A tibble: 3 × 9
## .y. group1 group2 n1 n2 statistic p p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <chr>
## 1 Liver_STm_CFUg GF E. clost… 7 11 45.5 0.554 0.807 ns
## 2 Liver_STm_CFUg GF A. faecis 7 5 15.5 0.807 0.807 ns
## 3 Liver_STm_CFUg GF E. coli 7 6 36.5 0.024 0.073 ns
vivo_cfu$Spleen_STm_CFUg <- as.numeric(vivo_cfu$Spleen_STm_CFUg)
## Warning: NAs introduced by coercion
vivo_cfu$Spleen_STm_CFUg[vivo_cfu$Spleen_STm_CFUg == 0] <- 1
p_spleen <- ggboxplot(vivo_cfu, x = "Group", y="Spleen_STm_CFUg",
fill = "Group",
palette = monocol_pal,
size = 1,
ylab = "STm CFU/g spleen",
xlab = "Spleen",
add = c("dotplot")) +
stat_compare_means(comparisons = monocol_comp,
label = "p.signif",
bracket.size = 0.5,
tip.length = 0)
p_spleen <- ggpar(p_spleen, legend = "right") + yscale("log10") +
scale_y_log10() +
theme(axis.ticks.x = element_blank(),
axis.text.x = element_blank(),
legend.title = element_blank())
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
p_spleen
## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_boxplot()`).
## Bin width defaults to 1/30 of the range of the data. Pick better value with
## `binwidth`.
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`stat_bindot()`).
## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_signif()`).
## Warning in wilcox.test.default(c(3.67669360962487, 3.54157924394658, 0, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(c(3.67669360962487, 3.54157924394658, 0, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(c(3.67669360962487, 3.54157924394658, 0, :
## cannot compute exact p-value with ties
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
# correct for multiple comparisons
wilcox_test(vivo_cfu, formula = Spleen_STm_CFUg ~ Group, comparisons = monocol_comp, p.adjust.method = "BH", ref.group = "GF")
## # A tibble: 3 × 9
## .y. group1 group2 n1 n2 statistic p p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <chr>
## 1 Spleen_STm_CFUg GF E. clos… 7 13 65 0.015 0.046 *
## 2 Spleen_STm_CFUg GF A. faec… 7 7 32 0.263 0.263 ns
## 3 Spleen_STm_CFUg GF E. coli 7 7 35 0.075 0.113 ns
pl2_p4 <- ggarrange(p_caecal_contents, p_mLN, p_liver, p_spleen, nrow = 1, align = "h",
common.legend = TRUE, legend = "right")
## Bin width defaults to 1/30 of the range of the data. Pick better value with
## `binwidth`.
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
## Bin width defaults to 1/30 of the range of the data. Pick better value with
## `binwidth`.
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
## Bin width defaults to 1/30 of the range of the data. Pick better value with
## `binwidth`.
## Warning in wilcox.test.default(c(1.07918124604762, 0, 0, 5.20411998265593, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(c(1.07918124604762, 0, 0, 5.20411998265593, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(c(1.07918124604762, 0, 0, 5.20411998265593, :
## cannot compute exact p-value with ties
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
## Warning: Removed 6 rows containing non-finite outside the scale range
## (`stat_boxplot()`).
## Bin width defaults to 1/30 of the range of the data. Pick better value with
## `binwidth`.
## Warning: Removed 6 rows containing missing values or values outside the scale range
## (`stat_bindot()`).
## Warning: Removed 6 rows containing non-finite outside the scale range
## (`stat_signif()`).
## Warning in wilcox.test.default(c(4.22788670461367, 1.2380461031288, 0,
## 2.46389298898591, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(c(4.22788670461367, 1.2380461031288, 0,
## 2.46389298898591, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(c(4.22788670461367, 1.2380461031288, 0,
## 2.46389298898591, : cannot compute exact p-value with ties
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_boxplot()`).
## Bin width defaults to 1/30 of the range of the data. Pick better value with
## `binwidth`.
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`stat_bindot()`).
## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_signif()`).
## Warning in wilcox.test.default(c(3.67669360962487, 3.54157924394658, 0, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(c(3.67669360962487, 3.54157924394658, 0, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(c(3.67669360962487, 3.54157924394658, 0, :
## cannot compute exact p-value with ties
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
pl2_p4
Colonisation resistance
# read data
cr <- read_excel(dp, "Colonisation_resistance")
# define names
cr$Group <- gsub(cr$Group, pattern = "A7", replacement = "E. coli") %>%
gsub(pattern = "A92", replacement = "E. clostridioformis") %>%
gsub(pattern = "A14", replacement = "A. faecis") %>%
gsub(pattern = "STm", replacement = "S. Typhimurium") %>%
gsub(pattern = "AnO2 S. Typhimurium", replacement = "S. Typhimurium")
# set levels
cr$Group <- factor(cr$Group, levels = unique(cr$Group)[c(3,5,2,4,1,7,6)])
# format data
cr$Bacterial_titres_CFU_g.mL <- as.numeric(cr$Bacterial_titres_CFU_g.mL)
cr$log_titres <- log10(cr$Bacterial_titres_CFU_g.mL)
# define colour palete
media_pal <- get_palette("Set1", 5)
names(media_pal) <- unique(cr$Group[cr$Assay == "Spent_STm"])[c(1,3,2,4,5)]
pspent <- ggline(cr[cr$Assay == "Spent_STm",], x = "Timepoint", y = "log_titres",
color = "Group",
palette = media_pal,
add = "mean_se",
numeric.x.axis = TRUE,
ylab = "S.Tm (CFU/mL)",
xlab = "Time (hrs)") +
scale_y_continuous(breaks = c(6,7,8,9), #limits = c(6,9),
labels = c("1e+06", "1e+07", "1e+08", "1e+09")) +
theme(legend.position = "right", legend.title = element_blank())
## Warning in stats::qt(ci/2 + 0.5, data_sum$length - 1): NaNs produced
pspent
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's fill values.
pstat <- ggline(cr[cr$Assay == "Stationary_STm",], x = "Timepoint", y = "log_titres",
color = "Group",
palette = media_pal,
add = c("mean_se"),
remove = 72,
numeric.x.axis = TRUE,
ylab = "S.Tm (CFU/mL)",
xlab = "Time (hrs)") +
scale_y_continuous(breaks = c(6,7,8,9),
labels = c("1e+06", "1e+07", "1e+08", "1e+09")) +
theme(legend.position = "right", legend.title = element_blank())
pstat
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's fill values.
p_comp <- ggarrange(pspent, pstat, nrow = 1, align = "hv",
common.legend = TRUE, legend = "right")
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's fill values.
## No shared levels found between `names(values)` of the manual scale and the
## data's fill values.
## No shared levels found between `names(values)` of the manual scale and the
## data's fill values.
p_comp
Commensal broth titres
com_titres <- cr[cr$Assay == "GrowthCurve_commensals" &
cr$Timepoint %in% c(24,36) &
cr$Group != "O2 S. Typhimurium",]
# summary statistics for titres
ddply(.data = com_titres, .variables = c("Group"),
.fun = function(x) {
c(summary(x$Bacterial_titres_CFU_g.mL),
mean_sd=sd(x$Bacterial_titres_CFU_g.mL))
})
## Group Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1 E. clostridioformis 8.50e+07 116750000 160000000 293625000 383000000 8.13e+08
## 2 A. faecis 2.85e+07 70875000 84850000 91837500 101250000 1.92e+08
## 3 E. coli 3.10e+08 480000000 682500000 656875000 786000000 1.08e+09
## 4 S. Typhimurium 2.45e+08 325000000 565000000 538571429 725000000 8.60e+08
## mean_sd
## 1 267401217
## 2 46935486
## 3 267959185
## 4 242344010
# significance of relationship
a_com <- aov(Bacterial_titres_CFU_g.mL ~ Group, com_titres)
TukeyHSD(a_com)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Bacterial_titres_CFU_g.mL ~ Group, data = com_titres)
##
## $Group
## diff lwr upr p adj
## A. faecis-E. clostridioformis -201787500 -510108705 106533705 0.2993976
## E. coli-E. clostridioformis 363250000 54928795 671571205 0.0163486
## S. Typhimurium-E. clostridioformis 244946429 -74196339 564089196 0.1785807
## E. coli-A. faecis 565037500 256716295 873358705 0.0001632
## S. Typhimurium-A. faecis 446733929 127591161 765876696 0.0036395
## S. Typhimurium-E. coli -118303571 -437446339 200839196 0.7424915
# plot data
p_ct <- ggbarplot(com_titres, x = "Group", y = "log_titres",
fill = "Group",
palette = media_pal,
add = c("mean_se"),
xlab = FALSE,
ylab = "Titres (CFU/mL)") +
scale_y_continuous(breaks = c(0,2,4,6,8,10), limits = c(0,10.2),
labels = c("1e+00", "1e+02", "1e+04", "1e+06", "1e+08", "1e+10")) +
theme(legend.position = "none") +
rotate_x_text(50)
p_ct
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
Supplementation with glucose
# get data
gluc <- read_excel(dp, "Spent_media_glucose")
# define names
gluc$Group <- gsub(gluc$Group, pattern = "A7", replacement = "E. coli") %>%
gsub(pattern = "A92", replacement = "E. clostridioformis") %>%
gsub(pattern = "A14", replacement = "A. faecis") %>%
gsub(pattern = "STm", replacement = "S. Typhimurium")
# define data groups to analyse
gluc_p <- gluc[gluc$Group %in% c("E. coli", "S. Typhimurium"),]
colnames(gluc_p)[c(3)] <- "Glucose"
# plot
p_gluc <- ggline(gluc_p,
x = "Glucose", y = "OD",
numeric.x.axis = TRUE,
xlab = "[Glucose], %",
ylab = "OD600",
add = "mean",
color = "Group",
palette = media_pal
) +
theme(legend.title = element_blank(), legend.position = "right")
p_gluc
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's fill values.
ex vivo faecal spiking assay
# get data
ex_vivo <- cr[cr$Assay == "Ex_vivo_STm",]
# plot
p_caecum <- ggline(ex_vivo[ex_vivo$Tissue_Broth == "Caecum",],
x = "Timepoint", y = "log_titres",
color = "Group", palette = get_palette("Set1", 4)[-3],
add = c("mean_se"),
numeric.x.axis = TRUE,
ylab = "S.Tm (CFU/g)",
xlab = "Time (hrs)") +
scale_y_continuous(breaks = c(6,7,8,9), labels = c("1e+06", "1e+07", "1e+08", "1e+09"))
## Warning in stats::qt(ci/2 + 0.5, data_sum$length - 1): NaNs produced
## Warning in stats::qt(ci/2 + 0.5, data_sum$length - 1): NaNs produced
## Warning in stats::qt(ci/2 + 0.5, data_sum$length - 1): NaNs produced
ggpar(p_caecum, legend = "none")
# significance calculations
a <- aov(Bacterial_titres_CFU_g.mL ~ Group, ex_vivo[ex_vivo$Assay == "Ex_vivo_STm" & ex_vivo$Tissue_Broth == "Caecum" & ex_vivo$Timepoint != 0,])
TukeyHSD(a)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Bacterial_titres_CFU_g.mL ~ Group, data = ex_vivo[ex_vivo$Assay == "Ex_vivo_STm" & ex_vivo$Tissue_Broth == "Caecum" & ex_vivo$Timepoint != 0, ])
##
## $Group
## diff lwr upr p adj
## E. clostridioformis-PBS -464666667 -1221102663 291769330 0.2986400
## E. coli-PBS -927981000 -1684416996 -171545004 0.0136425
## E. coli-E. clostridioformis -463314333 -900042859 -26585807 0.0358037
# read data
muc <- read_excel(path = dp, sheet = "Pre_STm_mucus_scores")
# curate data
muc$Expt <- lapply(strsplit(muc$Sample_index, split = "_"), function(x) {x[1]}) %>% unlist
muc$Mouse <- lapply(strsplit(muc$Sample_index, split = "_"), function(x) {x[2]}) %>% unlist
muc$Group <- NA
# define metadata
muc$Group <- gsub(pattern = "GP1 m1", replacement = "E. coli", paste(muc$Expt, muc$Mouse)) %>%
gsub(pattern = "GP1 M4", replacement = "E. coli") %>%
gsub(pattern = "GP1 M6", replacement = "E. clostridioformis") %>%
gsub(pattern = "GP1 M8", replacement = "E. clostridioformis") %>%
gsub(pattern = "GP2 M10", replacement = "E. clostridioformis") %>%
gsub(pattern = "GP2 M5", replacement = "E. coli") %>%
gsub(pattern = "GP3 M7", replacement = "E. coli") %>%
gsub(pattern = "GP3 M10", replacement = "E. clostridioformis") %>%
gsub(pattern = "GP2 M11", replacement = "E. clostridioformis") %>%
gsub(pattern = "GP3 M10", replacement = "E. clostridioformis") %>%
gsub(pattern = "DSS9 M1", replacement = "A. faecis") %>%
gsub(pattern = "DSS9 M3", replacement = "A. faecis") %>%
gsub(pattern = "DSS9 M2", replacement = "A. faecis") %>%
gsub(pattern = "GP3 M9", replacement = "E. clostridioformis") %>%
gsub(pattern = "GP2 M6", replacement = "E. coli") %>%
gsub(pattern = "DSS9 M4", replacement = "A. faecis") %>%
gsub(pattern = "GP3 M6", replacement = "E. coli")
muc$Group[grep("SPF", muc$Expt)] <- "SPF"
# define levels
muc$Group <- factor(muc$Group, levels = unique(muc$Group)[c(3,1,4,2)])
# caecum
# stats
wilcox_test(data = muc[muc$Tissue == "Caecum",], formula = Score ~ Group, ref.group = "SPF", p.adjust.method = "BH")
## # A tibble: 3 × 9
## .y. group1 group2 n1 n2 statistic p p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <chr>
## 1 Score SPF E. clostridio… 42 47 357 1.34e-7 4.02e-7 ****
## 2 Score SPF A. faecis 42 35 667 4.71e-1 4.71e-1 ns
## 3 Score SPF E. coli 42 42 653 3.4 e-2 5 e-2 ns
# plot
p <- ggviolin(muc[muc$Tissue == "Caecum",], x = "Group", y = "Score",
size = 1,
fill = "Group",
ylab = "Mucus penetration score",
palette = c(media_pal, "SPF" = "#A9A9A9")) +
scale_y_continuous(limits = c(0,5), breaks = c(0,1,2,3,4,5)) +
stat_compare_means(comparisons = list(c("SPF", "E. clostridioformis"),
c("E. clostridioformis", "A. faecis"),
c("E. clostridioformis", "E. coli")),
label = "p.signif", method = "wilcox", tip.length = 0) +
geom_boxplot(muc[muc$Tissue == "Caecum",], mapping = aes(x=Group, y=Score, fill=Group),
width = 0.2,
outliers = FALSE)
ggpar(p, xlab = FALSE, legend = "none") +
theme(legend.title = element_blank())
## Warning in wilcox.test.default(c(0, 2, 3, 0, 0, 1, 3, 0, 1, 0, 0, 1, 2, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(c(4, 3, 5, 5, 4, 1, 2, 3, 4, 5, 3, 1, 2, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(c(4, 3, 5, 5, 4, 1, 2, 3, 4, 5, 3, 1, 2, :
## cannot compute exact p-value with ties
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
## Warning: Removed 514 rows containing missing values or values outside the scale range
## (`geom_violin()`).
## Warning: Removed 9 rows containing missing values or values outside the scale range
## (`geom_signif()`).
# colon
# stats
wilcox_test(data = muc[muc$Tissue == "Colon",], formula = Score ~ Group, ref.group = "SPF", p.adjust.method = "BH")
## # A tibble: 3 × 9
## .y. group1 group2 n1 n2 statistic p p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <chr>
## 1 Score SPF E. clostridio… 36 38 177 4.97e-9 1.49e-8 ****
## 2 Score SPF A. faecis 36 14 42 1.9 e-7 1.9 e-7 ****
## 3 Score SPF E. coli 36 37 200. 2.95e-8 4.42e-8 ****
# plot
p <- ggviolin(muc[muc$Tissue == "Colon",], x = "Group", y = "Score",
size = 1,
fill = "Group",
#add = c("boxplot"), # manually specify below
ylab = "Mucus penetration score",
#title = "Caecum",
#add.params = list(fill = "lightgrey"),
palette = c(media_pal, "SPF" = "#A9A9A9")) +
scale_y_continuous(limits = c(0,5), breaks = c(0,1,2,3,4,5)) +
geom_boxplot(muc[muc$Tissue == "Colon",], mapping = aes(x=Group, y=Score, fill=Group),
width = 0.2,
outliers = FALSE)
ggpar(p, xlab = FALSE, legend = "none") +
theme(legend.title = element_blank())
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
## Warning: Removed 555 rows containing missing values or values outside the scale range
## (`geom_violin()`).
RTqPCR
d <- read_excel(path = dp, sheet = "IEC_RTqPCR_data")
d$Group <- gsub(d$Group, pattern = "A7", replacement = "E. coli") %>%
gsub(pattern = "A92", replacement = "E. clostridioformis") %>%
gsub(pattern = "A14", replacement = "A. faecis")
d$Group <- factor(d$Group, levels = unique(d$Group))
Barrier function target list: relmb reg3g saa1 cryptdin4 fut2 muc2 mZO-1 occludin nlrc4 il-1b
Immune function targets: ido1 il-6 raldh1 tnfa rdh7 tph1
# barrier function genes
barrier_genes <- c("relmb","reg3g","saa1","cryptdin4","fut2","muc2","mZO-1","occludin","nlrc4","il-1b")
d_bg <- d[d$Target %in% barrier_genes & d$Group != "S. muris" & d$Tissue %in% c("Caecum"),]
# establish levels
d_bg$Group <- factor(x = d_bg$Group, levels = unique(d_bg$Group)[c(1,2,5,3,4)])
# summary plots
lapply(split(d_bg, d_bg$Tissue), function(x) {
lapply(split(x, x$Target), function(y) {
p <- ggbarplot(y, x = "Group", y = "Expression",
add = c("mean_se"),
xlab = FALSE, ylab = "Relative Expression",
fill = "Group",
palette = c(monocol_pal, "SPF"="#A9A9A9"),
title = paste0(unique(y$Tissue), ": ", unique(y$Target)),) +
stat_compare_means(comparisons = list(c("GF", "E. coli"),
c("GF", "A. faecis"),
c("E. coli", "SPF"),
c("GF", "E. clostridioformis"),
c("A. faecis", "SPF"),
c("E. clostridioformis", "SPF"),
c("GF", "SPF")
),
label = "p.signif", tip.length = 0)
ggpar(p, legend = "right", x.text.angle = 30)
})
})
## $Caecum
## $Caecum$fut2
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
##
## $Caecum$`il-1b`
## Warning: Computation failed in `stat_signif()`.
## No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
## Caused by error in `wilcox.test.default()`:
## ! not enough 'y' observations
##
## $Caecum$muc2
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
##
## $Caecum$`mZO-1`
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
##
## $Caecum$nlrc4
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
##
## $Caecum$occludin
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
##
## $Caecum$reg3g
## Warning in wilcox.test.default(c(348.0103236912, 880.440031699341, 116.054073416104, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(c(348.0103236912, 880.440031699341, 116.054073416104, : No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
##
## $Caecum$relmb
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
##
## $Caecum$saa1
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
Caecum: fut2, muc2, mZO-1, nlrc4, occludin, reg3g, relmbm, saa1 - complete il-1b - no SPF -> exclude il-1b
# generate multiple comparisons
lapply(split(d_bg, d_bg$Tissue), function(x) {
lapply(split(x, x$Target), function(y) {
res.out <- t_test(formula = Expression ~ Group, data = x, p.adjust.method = "BH", ref.group = "GF")
res.out$Tissue <- unique(y$Tissue)
res.out$Target <- unique(y$Target)
return(res.out)
})
})
## $Caecum
## $Caecum$fut2
## # A tibble: 4 × 12
## .y. group1 group2 n1 n2 statistic df p p.adj p.adj.signif
## <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 Expression GF E. clo… 124 96 -0.0318 203. 0.975 0.975 ns
## 2 Expression GF A. fae… 124 36 -0.198 59.5 0.844 0.975 ns
## 3 Expression GF E. coli 124 99 -0.364 186. 0.716 0.975 ns
## 4 Expression GF SPF 124 73 -3.21 72.6 0.002 0.008 **
## # ℹ 2 more variables: Tissue <chr>, Target <chr>
##
## $Caecum$`il-1b`
## # A tibble: 4 × 12
## .y. group1 group2 n1 n2 statistic df p p.adj p.adj.signif
## <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 Expression GF E. clo… 124 96 -0.0318 203. 0.975 0.975 ns
## 2 Expression GF A. fae… 124 36 -0.198 59.5 0.844 0.975 ns
## 3 Expression GF E. coli 124 99 -0.364 186. 0.716 0.975 ns
## 4 Expression GF SPF 124 73 -3.21 72.6 0.002 0.008 **
## # ℹ 2 more variables: Tissue <chr>, Target <chr>
##
## $Caecum$muc2
## # A tibble: 4 × 12
## .y. group1 group2 n1 n2 statistic df p p.adj p.adj.signif
## <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 Expression GF E. clo… 124 96 -0.0318 203. 0.975 0.975 ns
## 2 Expression GF A. fae… 124 36 -0.198 59.5 0.844 0.975 ns
## 3 Expression GF E. coli 124 99 -0.364 186. 0.716 0.975 ns
## 4 Expression GF SPF 124 73 -3.21 72.6 0.002 0.008 **
## # ℹ 2 more variables: Tissue <chr>, Target <chr>
##
## $Caecum$`mZO-1`
## # A tibble: 4 × 12
## .y. group1 group2 n1 n2 statistic df p p.adj p.adj.signif
## <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 Expression GF E. clo… 124 96 -0.0318 203. 0.975 0.975 ns
## 2 Expression GF A. fae… 124 36 -0.198 59.5 0.844 0.975 ns
## 3 Expression GF E. coli 124 99 -0.364 186. 0.716 0.975 ns
## 4 Expression GF SPF 124 73 -3.21 72.6 0.002 0.008 **
## # ℹ 2 more variables: Tissue <chr>, Target <chr>
##
## $Caecum$nlrc4
## # A tibble: 4 × 12
## .y. group1 group2 n1 n2 statistic df p p.adj p.adj.signif
## <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 Expression GF E. clo… 124 96 -0.0318 203. 0.975 0.975 ns
## 2 Expression GF A. fae… 124 36 -0.198 59.5 0.844 0.975 ns
## 3 Expression GF E. coli 124 99 -0.364 186. 0.716 0.975 ns
## 4 Expression GF SPF 124 73 -3.21 72.6 0.002 0.008 **
## # ℹ 2 more variables: Tissue <chr>, Target <chr>
##
## $Caecum$occludin
## # A tibble: 4 × 12
## .y. group1 group2 n1 n2 statistic df p p.adj p.adj.signif
## <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 Expression GF E. clo… 124 96 -0.0318 203. 0.975 0.975 ns
## 2 Expression GF A. fae… 124 36 -0.198 59.5 0.844 0.975 ns
## 3 Expression GF E. coli 124 99 -0.364 186. 0.716 0.975 ns
## 4 Expression GF SPF 124 73 -3.21 72.6 0.002 0.008 **
## # ℹ 2 more variables: Tissue <chr>, Target <chr>
##
## $Caecum$reg3g
## # A tibble: 4 × 12
## .y. group1 group2 n1 n2 statistic df p p.adj p.adj.signif
## <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 Expression GF E. clo… 124 96 -0.0318 203. 0.975 0.975 ns
## 2 Expression GF A. fae… 124 36 -0.198 59.5 0.844 0.975 ns
## 3 Expression GF E. coli 124 99 -0.364 186. 0.716 0.975 ns
## 4 Expression GF SPF 124 73 -3.21 72.6 0.002 0.008 **
## # ℹ 2 more variables: Tissue <chr>, Target <chr>
##
## $Caecum$relmb
## # A tibble: 4 × 12
## .y. group1 group2 n1 n2 statistic df p p.adj p.adj.signif
## <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 Expression GF E. clo… 124 96 -0.0318 203. 0.975 0.975 ns
## 2 Expression GF A. fae… 124 36 -0.198 59.5 0.844 0.975 ns
## 3 Expression GF E. coli 124 99 -0.364 186. 0.716 0.975 ns
## 4 Expression GF SPF 124 73 -3.21 72.6 0.002 0.008 **
## # ℹ 2 more variables: Tissue <chr>, Target <chr>
##
## $Caecum$saa1
## # A tibble: 4 × 12
## .y. group1 group2 n1 n2 statistic df p p.adj p.adj.signif
## <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 Expression GF E. clo… 124 96 -0.0318 203. 0.975 0.975 ns
## 2 Expression GF A. fae… 124 36 -0.198 59.5 0.844 0.975 ns
## 3 Expression GF E. coli 124 99 -0.364 186. 0.716 0.975 ns
## 4 Expression GF SPF 124 73 -3.21 72.6 0.002 0.008 **
## # ℹ 2 more variables: Tissue <chr>, Target <chr>
# REG3G
t_test(data = d_bg[d_bg$Tissue == "Caecum" & d_bg$Target == "reg3g",], formula = Expression ~ Group, ref.group = "GF", p.adjust.method = "BH")
## # A tibble: 4 × 10
## .y. group1 group2 n1 n2 statistic df p p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 Expressi… GF E. cl… 14 11 1.24 21.2 2.28e-1 0.304 ns
## 2 Expressi… GF A. fa… 14 4 -4.54 10.4 9.76e-4 0.004 **
## 3 Expressi… GF E. co… 14 11 -0.582 10.8 5.73e-1 0.573 ns
## 4 Expressi… GF SPF 14 8 -3.77 7.92 6 e-3 0.011 *
# Expression GF E. clostridioformis 14 11 1.2414271 21.195713 0.228000 0.304 ns
# Expression GF A. faecis 14 4 -4.5413132 10.378775 0.000976 0.004 **
# Expression GF E. coli 14 11 -0.5817792 10.754973 0.573000 0.573 ns
# Expression GF SPF 14 8 -3.7701815 7.923493 0.006000 0.011 *
# RELMB
t_test(data = d_bg[d_bg$Tissue == "Caecum" & d_bg$Target == "relmb",], formula = Expression ~ Group, ref.group = "GF", p.adjust.method = "BH")
## # A tibble: 4 × 10
## .y. group1 group2 n1 n2 statistic df p p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 Expression GF E. clo… 14 11 -3.64 10.4 0.004 0.009 **
## 2 Expression GF A. fae… 14 4 -0.464 7.30 0.656 0.875 ns
## 3 Expression GF E. coli 14 11 -0.0709 15.6 0.944 0.944 ns
## 4 Expression GF SPF 14 10 -3.92 9.06 0.003 0.009 **
# Expression GF E. clostridioformis 14 11 -3.63797370 10.398726 0.004 0.009 **
# Expression GF A. faecis 14 4 -0.46391192 7.304980 0.656 0.875 ns
# Expression GF E. coli 14 11 -0.07089648 15.608042 0.944 0.944 ns
# Expression GF SPF 14 10 -3.92361910 9.057568 0.003 0.009 **
# SAA1
t_test(data = d_bg[d_bg$Tissue == "Caecum" & d_bg$Target == "saa1",], formula = Expression ~ Group, ref.group = "GF", p.adjust.method = "BH")
## # A tibble: 4 × 10
## .y. group1 group2 n1 n2 statistic df p p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 Expression GF E. clo… 14 11 2.55 23.0 0.018 0.072 ns
## 2 Expression GF A. fae… 14 4 -0.777 13.5 0.451 0.503 ns
## 3 Expression GF E. coli 14 11 -0.987 18.3 0.337 0.503 ns
## 4 Expression GF SPF 14 10 -0.694 10.6 0.503 0.503 ns
# Expression GF E. clostridioformis 14 11 2.5491894 22.99982 0.018 0.072 ns
# Expression GF A. faecis 14 4 -0.7766063 13.51839 0.451 0.503 ns
# Expression GF E. coli 14 11 -0.9870059 18.25768 0.337 0.503 ns
# Expression GF SPF 14 10 -0.6935388 10.64449 0.503 0.503 ns
# Tjp1
t_test(data = d_bg[d_bg$Tissue == "Caecum" & d_bg$Target == "mZO-1",], formula = Expression ~ Group, ref.group = "GF", p.adjust.method = "BH")
## # A tibble: 4 × 10
## .y. group1 group2 n1 n2 statistic df p p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 Expression GF E. clo… 14 11 -0.0888 15.2 0.93 0.93 ns
## 2 Expression GF A. fae… 14 4 1.30 13.3 0.215 0.43 ns
## 3 Expression GF E. coli 14 11 0.168 22.7 0.868 0.93 ns
## 4 Expression GF SPF 14 9 -2.18 8.27 0.059 0.238 ns
# Expression GF E. clostridioformis 14 11 -0.08880256 15.16124 0.930 0.930 ns
# Expression GF A. faecis 14 4 1.30078286 13.25484 0.215 0.430 ns
# Expression GF E. coli 14 11 0.16785856 22.69868 0.868 0.930 ns
# Expression GF SPF 14 9 -2.18372650 8.26687 0.059 0.238 ns
# generate publication plots
plot_rtqpcr <- function(tissue, target, ylh) {
df <- d_bg[d_bg$Tissue == tissue & d_bg$Target == target,]
p <- ggbarplot(df, x = "Group", y = "Expression",
add = c("mean_se"),
xlab = FALSE, ylab = "Relative Expression",
fill = "Group",
palette = c(monocol_pal, "SPF"="#A9A9A9"),
title = paste0(unique(df$Tissue), ": ", unique(df$Target)),) +
stat_compare_means(comparisons = list(c("GF", "E. coli"),
c("GF", "A. faecis"),
c("E. coli", "SPF"),
c("GF", "E. clostridioformis"),
c("A. faecis", "SPF"),
c("E. clostridioformis", "SPF"),
c("GF", "SPF")
),
label.y = ylh,
label = "p.signif", tip.length = 0, method = "t.test")
ggpar(p, legend = "none", x.text.angle = 30)
}
get_ylh_auto <- function(start) {
interval=start/7.5
c(start, start+interval, start+interval, start+(2*interval), start+(2*interval), start+(3*interval), start+(4*interval))
}
# caecum
plot_rtqpcr(tissue = "Caecum", target = "relmb", ylh = get_ylh_auto(start= 107000))
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
plot_rtqpcr(tissue = "Caecum", target = "reg3g", ylh = get_ylh_auto(start= 2600))
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
plot_rtqpcr(tissue = "Caecum", target = "saa1", ylh = get_ylh_auto(start= 135000))
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
plot_rtqpcr(tissue = "Caecum", target = "fut2", ylh = get_ylh_auto(start= 49000))
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
plot_rtqpcr(tissue = "Caecum", target = "muc2", ylh = get_ylh_auto(start= 930000))
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
plot_rtqpcr(tissue = "Caecum", target = "mZO-1", ylh = get_ylh_auto(start= 27500))
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
plot_rtqpcr(tissue = "Caecum", target = "occludin", ylh = get_ylh_auto(start= 5100))
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
plot_rtqpcr(tissue = "Caecum", target = "nlrc4", ylh = get_ylh_auto(start= 5400))
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
# get data
deg_df <- read.delim("curated_degs_BC_uniprot.tsv")
# summary statistics
ddply(.data = deg_df, .variables = c("Group"),
.fun = function(x){
data.frame(sig_degs=count(x$FDR <= 0.05),
nonsig_degs=count(x$FDR > 0.05))
})
## Group sig_degs.x sig_degs.freq nonsig_degs.x nonsig_degs.freq
## 1 A7 FALSE 14971 FALSE 65
## 2 A7 TRUE 65 TRUE 14971
## 3 A92 FALSE 14673 FALSE 363
## 4 A92 TRUE 363 TRUE 14673
## 5 SPF FALSE 13612 FALSE 1424
## 6 SPF TRUE 1424 TRUE 13612
# total genes
# A92 = 363
# A7 = 65
# SPF = 1424
# all three groups
names_s123 <- which(deg_df[deg_df$FDR <= 0.05 & !is.na(deg_df$FDR),1] %>% table ==3) %>% names
length(names_s123)
## [1] 6
# = 6
# pairwise sharing, excluding genes shared by all groups
# A92 vs A7
count(deg_df[deg_df$FDR <= 0.05 & !is.na(deg_df$FDR) & deg_df$Group == "A7" & !deg_df$Gene %in% names_s123,1] %in% deg_df[deg_df$FDR <= 0.05 & !is.na(deg_df$FDR) & deg_df$Group == "A92" & !deg_df$Gene %in% names_s123,1])
## x freq
## 1 FALSE 55
## 2 TRUE 4
# 4
# A92 vs SPF
count(deg_df[deg_df$FDR <= 0.05 & !is.na(deg_df$FDR) & deg_df$Group == "SPF"& !deg_df$Gene %in% names_s123,1] %in% deg_df[deg_df$FDR <= 0.05 & !is.na(deg_df$FDR) & deg_df$Group == "A92" & !deg_df$Gene %in% names_s123,1])
## x freq
## 1 FALSE 1279
## 2 TRUE 139
# 139
# SPF vs A7
count(deg_df[deg_df$FDR <= 0.05 & !is.na(deg_df$FDR) & deg_df$Group == "A7" & !deg_df$Gene %in% names_s123,1] %in% deg_df[deg_df$FDR <= 0.05 & !is.na(deg_df$FDR) & deg_df$Group == "SPF" & !deg_df$Gene %in% names_s123,1])
## x freq
## 1 FALSE 49
## 2 TRUE 10
# 10
# review significantly differentially expressed genes (DEGs)
# A7
deg_df[deg_df$Group == "A7" & deg_df$FDR <= 0.05 & !is.na(deg_df$FDR),]
# A92
deg_df[deg_df$Group == "A92" & deg_df$FDR <= 0.05 & !is.na(deg_df$FDR),]
# SPF
deg_df[deg_df$Group == "SPF" & deg_df$FDR <= 0.05 & !is.na(deg_df$FDR),]
deg_df$minuslog_padj <- -log(deg_df$FDR, base = 10)
# define colours by significance
deg_df$sig_status <- "lightgrey"
deg_df$sig_status[deg_df$FDR <= 0.05 & deg_df$log2FoldChange >= 1.5 & !is.na(deg_df$FDR)] <- "blue"
deg_df$sig_status[deg_df$FDR <= 0.05 & deg_df$log2FoldChange <= -1.5 & !is.na(deg_df$FDR)] <- "red"
colpal <- unique(deg_df$sig_status)
names(colpal) <- unique(deg_df$sig_status)
# volcano plots with labels
# E. coli
tmp_vp <- deg_df[deg_df$Group == "A7" & !is.na(deg_df$FDR),]
# add labels
for_labs <- tmp_vp[abs(tmp_vp$log2FoldChange) >= 1.5 & tmp_vp$FDR <= 0.05,]
labs_select <- unique(c(for_labs$Gene[for_labs$minuslog_padj > -log(0.05, base = 10) + (7.5*sd(tmp_vp$minuslog_padj))],
for_labs$Gene[abs(for_labs$log2FoldChange) >= 1.5+(5*sd(tmp_vp$log2FoldChange))]))
# plot
ggscatter(tmp_vp,
x = "log2FoldChange", y = "minuslog_padj",
title = "E. coli vs GF",
color = "sig_status",
xlab = "log2(Fold Change)",
ylab = "-log10(p-value)",
label = "Gene",
label.select = labs_select,
repel = TRUE,
palette = colpal) +
theme(legend.position = "none") +
scale_x_continuous(limits = c(-26,26))
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's fill values.
# A92
tmp_vp <- deg_df[deg_df$Group == "A92" & !is.na(deg_df$FDR),]
# labels
for_labs <- tmp_vp[abs(tmp_vp$log2FoldChange) >= 1.5 & tmp_vp$FDR <= 0.05,]
labs_select <- unique(c(for_labs$Gene[for_labs$minuslog_padj > -log(0.05, base = 10) + (7.5*sd(tmp_vp$minuslog_padj))],
for_labs$Gene[abs(for_labs$log2FoldChange) >= 1.5+(5*sd(tmp_vp$log2FoldChange))]))
# plot
ggscatter(tmp_vp,
x = "log2FoldChange", y = "minuslog_padj",
title = "E. clostridioformis vs GF",
color = "sig_status",
xlab = "log2(Fold Change)",
ylab = "-log10(p-value)",
label = "Gene",
label.select = labs_select,
repel = TRUE,
palette = colpal) +
theme(legend.position = "none") +
scale_x_continuous(limits = c(-32,32))
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's fill values.
# load QuickGO annotations for cell cycle
cc_go <- read_excel(path = dp, sheet = "GO_cell_cycle")
cc_go <- cc_go[which(cc_go$`GENE PRODUCT ID` %in% deg_df$ID),] %>% unique
# get significant cell cycle associated genes
cc_go_A92 <- deg_df[deg_df$Group == "A92" & deg_df$FDR <= 0.05 & !is.na(deg_df$ID) & !is.na(deg_df$FDR) & deg_df$ID %in% cc_go$`GENE PRODUCT ID`,]
cc_go_d <- deg_df[deg_df$Gene %in% cc_go_A92$Gene,]
cc_go_dm <- dcast(cc_go_d, Group ~ Gene, value.var = "log2FoldChange")
row.names(cc_go_dm) <- c("E. coli", "E. clostridioformis", "SPF")
cc_go_dm <- as.matrix(cc_go_dm[,-1])
# generate heatmap
pheatmap(mat = cc_go_dm, cluster_rows = F, cluster_cols = F, angle_col = 90,
color = rev(get_palette("RdBu", 1024)[c(100:924)]),
breaks = seq(-max(c(cc_go_dm)), max(c(cc_go_dm)), length.out = 825))
# load QuickGO annotations for leukocyte activation
la_go <- read_excel(path = dp, sheet = "GO_leukocyte_activation")
la_go <- la_go[which(la_go$`GENE PRODUCT ID` %in% deg_df$ID),] %>% unique
# get significant leukocyte activation associated genes
la_go_A7 <- deg_df[deg_df$Group == "A7" & deg_df$FDR <= 0.05 & !is.na(deg_df$ID) & !is.na(deg_df$FDR) & deg_df$ID %in% la_go$`GENE PRODUCT ID`,]
la_go_d <- deg_df[deg_df$Gene %in% la_go_A7$Gene,]
la_go_dm <- dcast(la_go_d, Group ~ Gene, value.var = "log2FoldChange")
row.names(la_go_dm) <- c("E. coli", "E. clostridioformis", "SPF")
la_go_dm <- as.matrix(la_go_dm[,-1])
# generate heatmap
pheatmap(mat = la_go_dm, cluster_rows = F, cluster_cols = F, angle_col = 90,
color = rev(get_palette("RdBu", 1024)[c(100:924)]),
breaks = seq(-max(c(la_go_dm)), max(c(la_go_dm)), length.out = 825))
Immune phenotyping and flow cytometry
# load data
iel_data <- read_excel(path = dp, sheet = "Flow_data_IEL")
noniel_data <- read_excel(path = dp, sheet = "Flow_data_nonIEL")
# CIEL data plots
d <- iel_data[iel_data$Tissue == "Caecum",]
# normalise to CD3+ instead of Viable
d$`gd-T cells | CD3+` <- d$`gd-T cells | Viable` / d$`CD3+ | Viable`*100
d$`CD8aa T cells | CD3+` <- d$`CD8aa T cells | Viable` / d$`CD3+ | Viable`*100
d$`CD8aa gd-T cells | CD3+` <- d$`CD8aa gd-T cells | Viable` / d$`CD3+ | Viable`*100
# make palette and order
monocol_pal2 <- c("#E41A1C","#377EB8","#4DAF4A","#984EA3")
names(monocol_pal2) <- c("GF", "E. clos.", "A. faecis", "E. coli")
iso_order2 <- c("GF","E. clos.","E. coli")
my_comparisons2 <- list(c("E. clos.","GF"), c("E. coli","GF"))
# build plots
df <- data.frame(Metadata=d$Group2, Var_group=d$`gd-T cells | CD3+`)
wilcox_test(data = df, formula = Var_group ~ Metadata, ref.group = "GF", p.adjust.method = "BH")
## # A tibble: 2 × 9
## .y. group1 group2 n1 n2 statistic p p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <chr>
## 1 Var_group GF E. clos. 11 11 83 0.151 0.151 ns
## 2 Var_group GF E. coli 11 12 22 0.006 0.011 *
# Var_group GF E. clos. 11 11 83 0.151 0.151 ns
# 2 Var_group GF E. coli 11 12 22 0.006 0.011 *
p_gdtc_cd3_ciel <- ggbarplot(df, x = "Metadata", y = "Var_group",
fill = "Metadata",
palette = monocol_pal2,
order = iso_order2,
ylab = "TCRg+ (% of CD3+)",
xlab = FALSE,
add = c("median_mad"),
add.params = list(width = 0.3)) +
stat_compare_means(comparisons = my_comparisons2, label = "p.signif") +
theme(legend.position = "none") +
rotate_x_text(30) +
geom_dotplot(binaxis = "y", stackdir = "center", binpositions="all", binwidth = max(df$Var_group)/50)
df <- data.frame(Metadata=d$Group2, Var_group=d$`CD8aa gd-T cells | CD3+`)
wilcox_test(data = df, formula = Var_group ~ Metadata, ref.group = "GF", p.adjust.method = "BH")
## # A tibble: 2 × 9
## .y. group1 group2 n1 n2 statistic p p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <chr>
## 1 Var_group GF E. clos. 11 11 75 0.358 0.358 ns
## 2 Var_group GF E. coli 11 12 23 0.009 0.018 *
# Var_group GF E. clos. 11 11 75 0.358 0.358 ns
# 2 Var_group GF E. coli 11 12 23 0.009 0.018 *
p_cd8aa_gdtc_cd3_ciel <- ggbarplot(df, x = "Metadata", y = "Var_group",
fill = "Metadata",
palette = monocol_pal2,
order = iso_order2,
ylab = "CD8aa+ TCRg+ (% of CD3+)",
xlab = FALSE,
add = c("median_mad"),
add.params = list(width = 0.3)) +
stat_compare_means(comparisons = my_comparisons2, label = "p.signif") +
theme(legend.position = "none") +
rotate_x_text(30) +
geom_dotplot(binaxis = "y", stackdir = "center", binpositions="all", binwidth = max(df$Var_group)/50)
df <- data.frame(Metadata=d$Group2, Var_group=d$`CD8aa T cells | CD3+`)
wilcox_test(data = df, formula = Var_group ~ Metadata, ref.group = "GF", p.adjust.method = "BH")
## # A tibble: 2 × 9
## .y. group1 group2 n1 n2 statistic p p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <chr>
## 1 Var_group GF E. clos. 11 11 78 0.27 0.54 ns
## 2 Var_group GF E. coli 11 12 69 0.88 0.88 ns
# Var_group GF E. clos. 11 11 78 0.27 0.54 ns
# 2 Var_group GF E. coli 11 12 69 0.88 0.88 ns
p_cd8aa_tc_cd3_ciel <- ggbarplot(df, x = "Metadata", y = "Var_group",
fill = "Metadata",
palette = monocol_pal2,
order = iso_order2,
ylab = "CD8aa+ TCRb+ (% of CD3+)",
xlab = FALSE,
add = c("median_mad"),
add.params = list(width = 0.3)) +
stat_compare_means(comparisons = my_comparisons2, label = "p.signif") +
theme(legend.position = "none") +
rotate_x_text(30) +
geom_dotplot(binaxis = "y", stackdir = "center", binpositions="all", binwidth = max(df$Var_group)/50)
ggarrange(p_gdtc_cd3_ciel, p_cd8aa_gdtc_cd3_ciel, p_cd8aa_tc_cd3_ciel,
ncol=3, nrow=1, align = "hv")
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
## Warning in wilcox.test.default(c(2.48574686431015, 1.05027932960894,
## 1.04510451045105, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(c(1.6270337922403, 1.72571428571429,
## 2.04200700116686, : cannot compute exact p-value with ties
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
## No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
# CLPL data plots
d <- noniel_data[noniel_data$Tissue == "Colon",]
df <- data.frame(Metadata=d$Group2, Var_group=d$`CD4 T cells | T cells`)
wilcox_test(data = df, formula = Var_group ~ Metadata, ref.group = "GF", p.adjust.method = "BH")
## # A tibble: 2 × 9
## .y. group1 group2 n1 n2 statistic p p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <chr>
## 1 Var_group GF E. clos. 11 11 38 0.151 0.151 ns
## 2 Var_group GF E. coli 11 12 33.5 0.049 0.098 ns
# Var_group GF E. clos. 11 10 22 0.022000 0.022000 *
# 2 Var_group GF E. coli 11 12 4 0.000153 0.000306 ***
p_cd4_tc_clpl <- ggbarplot(df, x = "Metadata", y = "Var_group",
fill = "Metadata",
palette = monocol_pal2,
order = iso_order2,
ylab = "CD4+ (% of TCRb+)",
xlab = FALSE,
add = c("median_mad"),
add.params = list(width = 0.3)) +
stat_compare_means(comparisons = my_comparisons2, label = "p.signif") +
theme(legend.position = "none") +
rotate_x_text(30) +
geom_dotplot(binaxis = "y", stackdir = "center", binpositions="all", binwidth = max(df$Var_group, na.rm = TRUE)/50)
df <- data.frame(Metadata=d$Group2, Var_group=d$`Tregs | T cells`)
wilcox_test(data = df, formula = Var_group ~ Metadata, ref.group = "GF", p.adjust.method = "BH")
## # A tibble: 2 × 9
## .y. group1 group2 n1 n2 statistic p p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <chr>
## 1 Var_group GF E. clos. 11 11 29.5 0.045 0.056 ns
## 2 Var_group GF E. coli 11 12 34.5 0.056 0.056 ns
# Var_group GF E. clos. 11 10 18 0.010 0.020 *
# 2 Var_group GF E. coli 11 12 73 0.689 0.689 ns
p_tregs_tc_clpl <- ggbarplot(df, x = "Metadata", y = "Var_group",
fill = "Metadata",
palette = monocol_pal2,
order = iso_order2,
ylab = "FoxP3+ (% of TCRb+)",
xlab = FALSE,
add = c("median_mad"),
add.params = list(width = 0.3)) +
stat_compare_means(comparisons = my_comparisons2,
label = "p.signif") +
theme(legend.position = "none") +
rotate_x_text(30) +
geom_dotplot(binaxis = "y", stackdir = "center", binpositions="all", binwidth = max(df$Var_group, na.rm = TRUE)/50)
df <- data.frame(Metadata=d$Group2, Var_group=d$Treg_CD4_ratio)
wilcox_test(data = df, formula = Var_group ~ Metadata, ref.group = "GF", p.adjust.method = "BH")
## # A tibble: 2 × 9
## .y. group1 group2 n1 n2 statistic p p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <chr>
## 1 Var_group GF E. clos. 11 11 34 0.088 0.176 ns
## 2 Var_group GF E. coli 11 12 45 0.211 0.211 ns
# Var_group GF E. clos. 11 10 26 0.043 0.043 *
# 2 Var_group GF E. coli 11 12 100 0.037 0.043 *
p_tregs_cd4_clpl <- ggbarplot(df, x = "Metadata", y = "Var_group",
fill = "Metadata",
palette = monocol_pal2,
order = iso_order2,
ylab = "CD4 Treg:CD4 Teff ratio",
xlab = FALSE,
add = c("median_mad"),
add.params = list(width = 0.3)) +
stat_compare_means(comparisons = my_comparisons2,
label = "p.signif") +
theme(legend.position = "none") +
rotate_x_text(30) +
geom_dotplot(binaxis = "y", stackdir = "center", binpositions="all", binwidth = max(df$Var_group, na.rm = TRUE)/50)
ggarrange(p_cd4_tc_clpl, p_tregs_tc_clpl, p_tregs_cd4_clpl,
ncol=3, nrow=1, align = "hv")
## Warning in wilcox.test.default(c(78.8, 73.9, 68.9, 72.9, 80, 86, 76.1, 78.5, :
## cannot compute exact p-value with ties
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
## Warning in wilcox.test.default(c(17.8, 11.8, 8.11, 5.84, 12, 30.5, 19.4, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(c(13.5, 15.8, 12.6, 7.69, 24.2, 16.5, 5.58, :
## cannot compute exact p-value with ties
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
## No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
# CIEL data plots
d <- iel_data[iel_data$Tissue == "Caecum",]
df <- data.frame(Metadata=d$Group2, Var_group=d$Treg_CD4_ratio)
wilcox_test(data = df, formula = Var_group ~ Metadata, ref.group = "GF", p.adjust.method = "BH")
## # A tibble: 2 × 9
## .y. group1 group2 n1 n2 statistic p p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <chr>
## 1 Var_group GF E. clos. 11 11 9 0.000275 0.00055 ***
## 2 Var_group GF E. coli 11 12 56 0.566 0.566 ns
# Var_group GF E. clos. 11 11 9 0.000275 0.00055 ***
# 2 Var_group GF E. coli 11 12 56 0.566000 0.56600 ns
p_tregs_cd4_ciel <- ggbarplot(df, x = "Metadata", y = "Var_group",
fill = "Metadata",
palette = monocol_pal2,
order = iso_order2,
ylab = "CD4 Treg:CD4 Teff ratio",
xlab = FALSE,
add = c("median_mad")) +
stat_compare_means(comparisons = my_comparisons2, label = "p.signif") +
theme(legend.position = "none") +
rotate_x_text(30) +
geom_dotplot(binaxis = "y", stackdir = "center", binpositions="all", binwidth = max(df$Var_group, na.rm = TRUE)/50)
# LIEL data plots
d <- iel_data[iel_data$Tissue == "Colon",]
df <- data.frame(Metadata=d$Group2, Var_group=d$Treg_CD4_ratio)
wilcox_test(data = df, formula = Var_group ~ Metadata, ref.group = "GF", p.adjust.method = "BH")
## # A tibble: 2 × 9
## .y. group1 group2 n1 n2 statistic p p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <chr>
## 1 Var_group GF E. clos. 11 10 24 0.029 0.059 ns
## 2 Var_group GF E. coli 11 12 44 0.19 0.19 ns
# Var_group GF E. clos. 11 10 24 0.029 0.059 ns
# 2 Var_group GF E. coli 11 12 44 0.190 0.190 ns
p_tregs_cd4_liel <- ggbarplot(df, x = "Metadata", y = "Var_group",
fill = "Metadata",
palette = monocol_pal2,
order = iso_order2,
ylab = "CD4 Treg:CD4 Teff ratio",
xlab = FALSE,
add = c("median_mad")) +
stat_compare_means(comparisons = my_comparisons2, label = "p.signif") +
theme(legend.position = "none") +
rotate_x_text(30) +
geom_dotplot(binaxis = "y", stackdir = "center", binpositions="all", binwidth = max(df$Var_group, na.rm = TRUE)/50)
# SIEL data plots
d <- iel_data[iel_data$Tissue == "Small intestine",]
df <- data.frame(Metadata=d$Group2, Var_group=d$Treg_CD4_ratio)
wilcox_test(data = df, formula = Var_group ~ Metadata, ref.group = "GF", p.adjust.method = "BH")
## # A tibble: 2 × 9
## .y. group1 group2 n1 n2 statistic p p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <chr>
## 1 Var_group GF E. clos. 11 6 12 0.036 0.073 ns
## 2 Var_group GF E. coli 11 8 52 0.545 0.545 ns
# Var_group GF E. clos. 11 6 12 0.036 0.073 ns
# 2 Var_group GF E. coli 11 8 52 0.545 0.545 ns
p_tregs_cd4_siel <- ggbarplot(df, x = "Metadata", y = "Var_group",
fill = "Metadata",
palette = monocol_pal2,
order = iso_order2,
ylab = "CD4 Treg:CD4 Teff ratio",
xlab = FALSE,
add = c("median_mad")) +
stat_compare_means(comparisons = my_comparisons2, label = "p.signif") +
theme(legend.position = "none") +
rotate_x_text(30) +
geom_dotplot(binaxis = "y", stackdir = "center", binpositions="all", binwidth = max(df$Var_group, na.rm = TRUE)/50)
ggarrange(p_tregs_cd4_ciel, p_tregs_cd4_liel, p_tregs_cd4_siel,
ncol=3, nrow=1, align = "hv")
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
## No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
## No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
# load data
spf <- read_excel(path = dp, sheet = "SPF_screen_data")
gf <- read_excel(path = dp, sheet = "GF_screen_data")
# rename treatment groups
spf$Abx[spf$Treatment == "PBS-no_abx"] <- "SPF"
spf$Abx[spf$Abx == "OG AVMN+Cf"] <- "Abx i.g."
spf$Abx[spf$Abx == "DW AVMN"] <- "Abx d.w."
# curate data
spf <- spf[spf$Commensal_taxonomy_r95 == "Control",]
gf <- gf[gf$Treatment == "PBS" & gf$Col_status == "Clean",]
# SPF
spf_perc <- apply(X = spf[,grep(pattern = "Weight_", colnames(spf))],
MARGIN = 1,
function(x) {
x <- as.numeric(x)
x/x[1]*100
}) %>% t %>% as.data.frame
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
colnames(spf_perc) <- paste0("dpi",seq(0,length(grep(pattern = "Weight_", colnames(spf)))-1))
spf <- cbind(spf, spf_perc)
# GF
gf_perc <- apply(X = gf[,grep(pattern = "Weight_", colnames(gf))],
MARGIN = 1,
function(x) {
x <- as.numeric(x)
x/x[1]*100
}) %>% t %>% as.data.frame
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
colnames(gf_perc) <- paste0("dpi",seq(0,length(grep(pattern = "Weight_", colnames(gf)))-1))
gf <- cbind(gf, gf_perc)
# setfill blanks to enable copmarison
gf$dpi5 <- NA
gf$dpi6 <- NA
gf$dpi7 <- NA
gf$dpi8 <- NA
gf$dpi9 <- NA
gf$dpi10 <- NA
gf$Abx <- "GF"
gf$Commensal_faeces_post_abx_d0 <- NA
gf$Commensal_faeces_post_abx_d1 <- NA
gf$Commensal_faeces_post_abx_d2 <- NA
# combine datasets
spf_gf <- rbind(gf[,colnames(gf) %in% colnames(spf)], spf[,colnames(spf) %in% colnames(gf)])
# set levels
spf_gf$Abx <- factor(spf_gf$Abx, levels = unique(spf_gf$Abx)[c(1,3,2,4)])
# prepare data
spf_gf$idvar <- paste0("M", seq(1, nrow(spf_gf), 1))
spf_gfm <- melt(spf_gf, measure.vars = colnames(spf_gf)[grep(pattern = "^dpi", colnames(spf_gf))])
spf_gfm$dpi <- gsub(pattern = "dpi", replacement = "", spf_gfm$variable) %>% as.numeric
# remove data when mice in a group has died
spf_gfm_v <- ddply(.data = spf_gfm, .variables = c("Abx", "dpi"),
.fun = function(x, y) {
i <- all(!is.na(x$value))
if (i == TRUE) {
dtmp <- x
} else {
dtmp <- NULL
}
return(dtmp)
})
# plot only up to day 3
p <- ggline(spf_gfm_v, x = "dpi", y = "value",
group = "Abx",
color = "Abx",
size = 1,
numeric.x.axis = TRUE,
plot_type = "b",
add = "mean_se",
palette = "Dark2",
xlab = "Days post-infection",
ylab = "Weight (% of day 0)") +
geom_hline(yintercept = 80, linetype = "dashed") +
theme(legend.title = element_blank()) +
scale_x_continuous(breaks = c(0,2,4,6,8))
ggpar(p, legend = "right")
# anova across timeoints
summary(a <- aov(value ~ Abx, spf_gfm_v))
## Df Sum Sq Mean Sq F value Pr(>F)
## Abx 3 1683 560.9 15.3 1.68e-08 ***
## Residuals 121 4437 36.7
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(a)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = value ~ Abx, data = spf_gfm_v)
##
## $Abx
## diff lwr upr p adj
## Abx d.w.-GF 2.822315 -1.731480 7.376110 0.3743277
## Abx i.g.-GF 5.261503 1.716378 8.806627 0.0010180
## SPF-GF 10.201133 6.185063 14.217202 0.0000000
## Abx i.g.-Abx d.w. 2.439188 -1.974457 6.852832 0.4772597
## SPF-Abx d.w. 7.378818 2.578697 12.178939 0.0006163
## SPF-Abx i.g. 4.939630 1.083204 8.796056 0.0061189
# wilcoxon test for weight loss at 1 dpi
spf_gfm_v[spf_gfm_v$dpi == 1,] %>%
wilcox_test(value ~ Abx) %>%
adjust_pvalue(method = "BH")
## # A tibble: 6 × 9
## .y. group1 group2 n1 n2 statistic p p.adj p.adj.signif
## <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <chr>
## 1 value GF Abx d.w. 18 9 29 6 e-3 9 e-3 *
## 2 value GF Abx i.g. 18 11 0 5.78e-8 3.47e-7 ****
## 3 value GF SPF 18 3 0 2 e-3 4 e-3 **
## 4 value Abx d.w. Abx i.g. 9 11 0 1.19e-5 3.57e-5 ****
## 5 value Abx d.w. SPF 9 3 0 9 e-3 1.08e-2 *
## 6 value Abx i.g. SPF 11 3 23 3.68e-1 3.68e-1 ns
## survival
spf_gf$Status <- 2
spf_gf$Status[spf_gf$Abx == "SPF"] <- 1
# format data
spf_gf$Survival_day <- apply(X=spf_gf[,grep(pattern = "^dpi", colnames(spf_gf))], MARGIN = 1, function(x) {
colnames(spf_gf)[grep(pattern = "^dpi", colnames(spf_gf))][length(x) + 1 - as.numeric(which(!is.na((rev(x))))[1])]
}) %>% gsub(pattern = "dpi", replacement = "") %>% as.numeric
spf_gf_survival <- ddply(.data = spf_gf, .variables = "Abx",
.fun = function(x) {
dpi0_surv <- length(which(!is.na(x$dpi1)))
dpi1_surv <- length(which(!is.na(x$dpi2)))
dpi2_surv <- length(which(!is.na(x$dpi3)))
dpi3_surv <- length(which(!is.na(x$dpi4)))
dpi4_surv <- length(which(!is.na(x$dpi5)))
dpi5_surv <- length(which(!is.na(x$dpi6)))
dpi6_surv <- length(which(!is.na(x$dpi7)))
dpi7_surv <- length(which(!is.na(x$dpi8)))
dpi8_surv <- length(which(!is.na(x$dpi9)))
dpi9_surv <- length(which(!is.na(x$dpi10)))
dpi0 <- dpi0_surv/dpi0_surv*100
dpi1 <- dpi1_surv/dpi0_surv*100
dpi2 <- dpi2_surv/dpi0_surv*100
dpi3 <- dpi3_surv/dpi0_surv*100
dpi4 <- dpi4_surv/dpi0_surv*100
dpi5 <- dpi5_surv/dpi0_surv*100
dpi6 <- dpi6_surv/dpi0_surv*100
dpi7 <- dpi7_surv/dpi0_surv*100
dpi8 <- dpi8_surv/dpi0_surv*100
dpi9 <- dpi9_surv/dpi0_surv*100
c(dpi0=dpi0,
dpi1=dpi1,
dpi2=dpi2,
dpi3=dpi3,
dpi4=dpi4,
dpi5=dpi5,
dpi6=dpi6,
dpi7=dpi7,
dpi8=dpi8,
dpi9=dpi9,
dpi10=0)
}
)
spf_gf_survival[,grep("dpi", colnames(spf_gf_survival))] <- apply(spf_gf_survival[,grep("dpi", colnames(spf_gf_survival))],
MARGIN = 1,
FUN = function(x) {
i <- sum(x == 0) - 1
if ( i != 0) {
x[12-c(1:i)] <- NA
}
as.numeric(x)
}) %>% t
spf_gfm_survival <- melt(data = spf_gf_survival, id.vars = "Abx")
spf_gfm_survival$dpi <- gsub(pattern = "dpi", replacement = "", spf_gfm_survival$variable) %>% as.numeric
spf_gfm_survival$value[spf_gfm_survival$Abx == "SPF" & spf_gfm_survival$variable == "dpi8"] <- 100
# plot only up to 8 days
spf_gfm_survival <- spf_gfm_survival[!spf_gfm_survival$variable %in% c("dpi9", "dpi10"),]
# plot survival
p <- ggplot(data = spf_gfm_survival, mapping = aes(x = dpi, y = value, group = Abx)) +
geom_step(aes(color = Abx), size = 1) +
#geom_point() +
theme_pubr() +
theme(legend.title = element_blank())
ggpar(p, palette = "Dark2", xlab = "Days post-infection", ylab = "Percent survival",
legend = "right") +
scale_x_continuous(breaks = c(0,2,4,6,8))
## Warning: Removed 10 rows containing missing values or values outside the scale range
## (`geom_step()`).
# load data
spf <- read_excel(path = dp, sheet = "SPF_screen_data")
spf <- spf[spf$Treatment != "PBS-no_abx",] # only keep antibiotic treated mice
# format data
spf_perc <- apply(X = spf[,grep(pattern = "Weight_", colnames(spf))],
MARGIN = 1,
function(x) {
x <- as.numeric(x)
x/x[1]*100
}) %>% t %>% as.data.frame
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
## Warning in FUN(newX[, i], ...): NAs introduced by coercion
colnames(spf_perc) <- paste0("dpi",seq(0,length(grep(pattern = "Weight_", colnames(spf)))-1))
spf <- cbind(spf, spf_perc)
spf$Commensal_faeces_post_abx_d0 <- as.numeric(spf$Commensal_faeces_post_abx_d0)
spf$Commensal_faeces_post_abx_d1 <- as.numeric(spf$Commensal_faeces_post_abx_d1)
## Warning: NAs introduced by coercion
spf$Commensal_faeces_post_abx_d2 <- as.numeric(spf$Commensal_faeces_post_abx_d2)
## Warning: NAs introduced by coercion
# calculate survival
spf$Survival_day <- apply(X=spf[,grep(pattern = "^dpi", colnames(spf))], MARGIN = 1, function(x) {
colnames(spf)[grep(pattern = "^dpi", colnames(spf))][length(x) + 1 - as.numeric(which(!is.na((rev(x))))[1])]
}) %>% gsub(pattern = "dpi", replacement = "") %>% as.numeric
# fix 0s post-log scale
spf$Log_d2_commensal <- log10(spf$Commensal_faeces_post_abx_d2) %>% gsub(pattern = "-Inf", replacement = 0) %>% as.numeric
# only include samples with correct data
spf2 <- spf[which(!is.na(spf$Commensal_faeces_post_abx_d2)),]
spf3 <- rbind(spf2[spf2$Treatment %in% c("PBS"),], # only use mice that were not treated with commensals prior to faeces plating
spf2[spf2$Experiment == "SPF_dw",]) %>% unique
p_comm_surv <- ggscatter(spf3, x = "Log_d2_commensal", y = "Survival_day",
add = "reg.line",
#add = "loess",
conf.int = TRUE, cor.coef = TRUE,
xlab = "CFU/g faeces, 2 days post-antibiotics",
ylab = "Survival (d.p.i.)") +
scale_x_continuous(breaks = c(0,2,4,6,8,10),
labels = c("N.D.", "1e+02", "1e+04", "1e+06", "1e+08", "1e+10"),
limits = c(0,10.5))
p_comm_surv
# format data
spf2m <- melt(spf2[spf2$Experiment == "SPF_dw",], measure.vars = colnames(spf2)[grep(colnames(spf2), pattern = "Commensal_faeces")][c(1:3)])
spf2m$dp_abx <- gsub(spf2m$variable, pattern = "Commensal_faeces_post_abx_d", replacement = "") %>% as.numeric
# correct 0s post log
spf2m$log_value <- log10(spf2m$value) %>% gsub(pattern = "-Inf", replacement = 0) %>% as.numeric
spf2m$Cage <- spf2m$Treatment
i=1
for (Cage in unique(spf2m$Treatment) %>% sort) {
spf2m$Cage[spf2m$Treatment == Cage] <- paste0("Cage", " ", i)
i=i+1
}
p <- ggbarplot(spf2m, x = "dp_abx", y = "log_value", group = "Treatment", add = "mean_se",
fill = "Cage",
position = position_dodge(0.9),
palette = "Greys",
xlab = "Days post-antibiotics",
ylab = "S.Tm CFU/g faeces") +
scale_y_continuous(breaks = c(0,2,4,6,8), labels = c("N.D.", "1e+02", "1e+04", "1e+06", "1e+08"))
## Warning in stats::qt(ci/2 + 0.5, data_sum$length - 1): NaNs produced
## Warning in stats::qt(ci/2 + 0.5, data_sum$length - 1): NaNs produced
## Warning in stats::qt(ci/2 + 0.5, data_sum$length - 1): NaNs produced
ggpar(p, legend = "right", legend.title = "")
# automate plot building with corrected p-values
build.IEL_plot <- function(i, site, ylab) {
d <- iel_data[iel_data$Tissue == site,]
i=i+4
print(colnames(d)[i])
df <- data.frame(d$Group2, d[,i])
group_name <- as.character(colnames(df)[2])
colnames(df) <- c("Metadata", "Var_group")
stat.test <- wilcox_test(data = df, formula = Var_group ~ Metadata, ref.group = "GF", p.adjust.method = "BH")
stat.test <- stat.test %>% mutate(y.position = c(max(df$Var_group, na.rm = TRUE)+(max(df$Var_group, na.rm = TRUE)/15),
max(df$Var_group, na.rm = TRUE)+(2*(max(df$Var_group, na.rm = TRUE)/15))))
p <- ggbarplot(df,
x = "Metadata",
y = "Var_group",
fill = "Metadata",
palette = monocol_pal2,
order = iso_order2,
ylab = ylab,
xlab = FALSE,
add = c("median_mad")) +
stat_pvalue_manual(stat.test, label = "p.adj") +
theme(legend.position = "none") +
rotate_x_text(30) +
geom_dotplot(binaxis = "y", stackdir = "center", binpositions="all", binwidth = max(df$Var_group, na.rm = TRUE)/50)
return(p)
}
build.nonIEL_plot <- function(i, site, ylab) {
d <- noniel_data[noniel_data$Tissue == site,]
i=i+4
print(colnames(d)[i])
df <- data.frame(d$Group2, d[,i])
group_name <- as.character(colnames(df)[2])
colnames(df) <- c("Metadata", "Var_group")
stat.test <- wilcox_test(data = df, formula = Var_group ~ Metadata, ref.group = "GF", p.adjust.method = "BH")
stat.test <- stat.test %>% mutate(y.position = c(max(df$Var_group, na.rm = TRUE)+(max(df$Var_group, na.rm = TRUE)/15),
max(df$Var_group, na.rm = TRUE)+(2*(max(df$Var_group, na.rm = TRUE)/15))))
p <- ggbarplot(df,
x = "Metadata",
y = "Var_group",
fill = "Metadata",
palette = monocol_pal2,
order = iso_order2,
ylab = ylab,
xlab = FALSE,
add = c("median_mad")) +
stat_pvalue_manual(stat.test, label = "p.adj") +
theme(legend.position = "none") +
rotate_x_text(30) +
geom_dotplot(binaxis = "y", stackdir = "center", binpositions="all", binwidth = max(df$Var_group, na.rm = TRUE)/50)
return(p)
}
# colon IEL
p_l1 <- build.IEL_plot(i = 3, site = "Colon", ylab = "B220+ (% of Viable CD45+)")
## [1] "B cells | Viable"
p_l2 <- build.IEL_plot(i = 15, site = "Colon", ylab = "CD8ab+ (% of TCRb+)")
## [1] "CD8ab T cells | T cells"
p_l3 <- build.IEL_plot(i = 14, site = "Colon", ylab = "CD4+ (% of TCRb+)")
## [1] "CD4 T cells | T cells"
p_l4 <- build.IEL_plot(i = 20, site = "Colon", ylab = "FoxP3+ (% of TCRb+)")
## [1] "Tregs | T cells"
# colon LPL
p_l5 <- build.nonIEL_plot(i = 3, site = "Colon", ylab = "B220+ (% of Viable CD45+)")
## [1] "B cells | Viable"
p_l6 <- build.nonIEL_plot(i = 10, site = "Colon", ylab = "CD8ab+ (% of TCRb+)")
## [1] "CD8ab T cells | T cells"
p_l7 <- build.nonIEL_plot(i = 9, site = "Colon", ylab = "CD4+ (% of TCRb+)")
## [1] "CD4 T cells | T cells"
p_l8 <- build.nonIEL_plot(i = 12, site = "Colon", ylab = "FoxP3+ (% of TCRb+)")
## [1] "Tregs | T cells"
# mLN
p_m1 <- build.nonIEL_plot(i = 3, site = "mLN", ylab = "B220+ (% of Viable CD45+)")
## [1] "B cells | Viable"
p_m2 <- build.nonIEL_plot(i = 10, site = "mLN", ylab = "CD8ab+ (% of TCRb+)")
## [1] "CD8ab T cells | T cells"
p_m3 <- build.nonIEL_plot(i = 9, site = "mLN", ylab = "CD4+ (% of TCRb+)")
## [1] "CD4 T cells | T cells"
p_m4 <- build.nonIEL_plot(i = 12, site = "mLN", ylab = "FoxP3+ (% of TCRb+)")
## [1] "Tregs | T cells"
ggarrange(p_l1, p_l5, p_m1,
p_l2, p_l6, p_m2,
p_l3, p_l7, p_m3,
p_l4, p_l8, p_m4,
nrow = 4, ncol = 3)
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
## No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
## No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
## No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
## No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
## No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
## No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
## No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
## No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
## No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
## No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
## No shared levels found between `names(values)` of the manual scale and the
## data's colour values.